function feature = colorJointPyramid(grayLayer, ALayer, BLayer, numBins,numLevels)
    pyramid = cell(2^(numLevels-1),2^(numLevels-1),numLevels);
    % we assume the images are of in gray level [0 1];.
    [sx sy]=size(grayLayer); % assume the three layers are of the same size;
    % for level = numLevels

    % for level numLevels
    for b_x = 1:2^(numLevels-1)
        for b_y = 1:2^(numLevels-1)
            x_block_size = ceil(sx*0.5^(numLevels-1));
            y_block_size = ceil(sy*0.5^(numLevels-1));
            x = x_block_size*(b_x-1)+1:min(sx,x_block_size*b_x);
            y = y_block_size*(b_y-1)+1:min(sy,y_block_size*b_y);
            counts=jointHistogram(grayLayer(x,y),ALayer(x,y),BLayer(x,y),numBins);
            pyramid{b_x,b_y,numLevels}=counts;
        end
    end

        
    for iLevel = numLevels-1:-1:1
        for bx = 1:2^(iLevel-1)
            for by = 1:2^(iLevel-1)
                counts = pyramid{bx*2-1,by*2-1,iLevel+1};
                counts = counts + pyramid{bx*2,by*2-1,iLevel+1};
                counts = counts + pyramid{bx*2-1,by*2,iLevel+1};
                counts = counts + pyramid{bx*2,by*2,iLevel+1};
                pyramid{bx,by,iLevel}=counts/4;
            end
        end
    end

    feature = [];
    for iLevel = numLevels:-1:1
        for b_x = 1:2^(iLevel-1)
            for b_y = 1:2^(iLevel-1)
                feature = [feature;pyramid{b_x,b_y,iLevel}];
            end
        end
    end


    % if have find scaled levels, will comput the histogram intersection between outer cells and inner cells of an image.
    if numLevels>=3
        outter_hist = zeros(prod(numBins),1);
        inner_hist = outter_hist;
        for j = 1:4
    	outter_hist = outter_hist + pyramid{1,j,3} + pyramid{4,j,3};
        end
        outter_hist = outter_hist + pyramid{2,1,3} + pyramid{2,4,3};
        outter_hist = outter_hist + pyramid{3,1,3} + pyramid{3,4,3};
        outter_hist = outter_hist/12;
        outter_hist = outter_hist(1:numBins(1));

        corner_hist = pyramid{1,1,3}+pyramid{4,4,3} + pyramid{1,4,3} + pyramid{4,1,3};  
        corner_hist = corner_hist(1:numBins(1))/4; 


        inner_hist = inner_hist + pyramid{2,2,3} + pyramid{2,3,3} + pyramid{3,2,3} + pyramid{3,3,3};
        inner_hist = inner_hist/4;
        inner_hist = inner_hist(1:numBins(1));

        feature = [feature; reshape(min(inner_hist,corner_hist),[],1)];
    end

end % function

function hist = jointHistogram(L,A,B,numBins)
    hist = zeros(prod(numBins),1);

    idxL = ceil(double(L+1e-6)/(1/numBins(1)+1e-6));
    idxA = ceil(double(A+1e-6)/(1/numBins(2)+1e-6));
    idxB = ceil(double(B+1e-6)/(1/numBins(3)+1e-6));

    idx = sub2ind(numBins,idxL,idxA,idxB);

    if length(idx)>prod(numBins)
        for iBin = 1:prod(numBins)
            hist(iBin)=length(find(idx==iBin));
        end
    else
        for i_pix = 1:length(idx)
            hist(idx(i_pix))  = hist(idx(i_pix)) + 1;
        end 
    end
    hist = hist/length(idx);
    hist = max(1e-3,hist);
end


